Plot the WASP14 spectrum with Regions overlaid
In [31]:
%matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter as FSF
from matplotlib.ticker import MaxNLocator
from matplotlib.ticker import MultipleLocator
import numpy as np
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed
In [7]:
from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import TRES, SPEX, HDF5Interface
import scipy.sparse as sp
myDataSpectrum = DataSpectrum.open("../../data/WASP14/WASP14-2009-06-14.hdf5", orders=np.array([22]))
myInstrument = TRES()
# myInstrument = SPEX()
myHDF5Interface = HDF5Interface("../../libraries/PHOENIX_TRES_F.hdf5")
#Load a model using the JSON file
myModel = Model.from_json("../WASP14_22_model_final_region.json", myDataSpectrum, myInstrument, myHDF5Interface)
In [3]:
myOrderModel = myModel.OrderModels[0]
model_flux = myOrderModel.get_spectrum()
In [8]:
spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
model_fl = myOrderModel.get_spectrum()
residuals = fl - model_fl
mask = spec.masks[0]
cov = myModel.OrderModels[0].get_Cov().todense()
In [ ]:
#also, we should visualize what typical draws look like from the Matern kernel
#For example, what does a random draw from Matern look like, what does 3 times that in amplitude look like...
In [6]:
N = len(wl)
fake_cov = np.random.multivariate_normal(np.zeros((N,)), cov)
In [ ]:
for cov in cov_list:
plt.plot(wl, cov, "k", lw=0.2)
plt.plot(wl, model_flux - fl, "b")
plt.show()
In [26]:
import StellarSpectra.constants as C
def gaussian(ax, region):
'''
Paint a Gaussian envelope on the axes.
'''
if type(region) == dict:
a = 10**region_dict['loga']
mu = region_dict['mu']
sigma = region_dict['sigma']
elif type(region) == list:
a, mu, sigma = region
a = 10**a
sig = sigma/C.c_kms * mu #convert to angstroms for the spacing
wl = np.linspace(mu - 4 * sig, mu + 4 * sig)
ys = a * np.exp(-0.5 * (C.c_kms/mu)**2 * (wl - mu)**2/sigma**2)
ax.fill_between(wl, -ys, ys, alpha=0.8, color="magenta", zorder=100)
In [15]:
#List of where the regions are. From 8/13/14 4sig run
region_list = [[ -1.37557276e+01, 5.17068109e+03, 5.07069537e+00],
[ -13.7420562, 5198.61903481, 5.48643032],
[ -1.34095891e+01, 5.23530766e+03, 4.82242126e+00],
[ -13.64347323, 5150.74527037, 6.67246923],
[ -13.70379713, 5169.07086766, 6.47508788],
[ -13.76186598, 5202.25227836, 7.45325824],
[ -13.70503714, 5187.83335919, 5.30911125],
[ -13.38718559, 5142.85381703, 7.06802524],
[ -13.58059031, 5151.82322042, 6.15930413],
[ -13.37354867, 5188.75488375, 5.79597804],
[ -13.74135283, 5234.57274185, 5.49249104]]
global_cov = 10**(-14.34)
In [42]:
fig, ax = plt.subplots(nrows=2, figsize=(6.5,3), sharex=True)
l0, = ax[0].plot(wl, fl, "b")
l1, = ax[0].plot(wl, model_flux, "r")
l2, = ax[1].plot(wl, residuals, color="0.5")
ax[0].legend([l0, l1], ["data", "model"], loc="lower right", ncol=2, prop={'size':10})
ax[1].set_xlabel(r"$\lambda$ [\AA]")
thresh = 3 * global_cov
ax[1].axhline(thresh, color="0.5", ls="-.")
ax[1].axhline(-thresh, color="0.5", ls="-.")
for region in region_list:
gaussian(ax[1], region)
ax[1].xaxis.set_major_formatter(FSF("%.0f"))
ax[1].set_xlim(np.min(wl), np.max(wl))
scale = 1.1 * np.max(np.abs(residuals))
ax[1].set_ylim(-scale, scale)
l3, = ax[1].plot(wl, residuals + 2, color="magenta") #Fake plot to just get the color handle for the legend
ax[1].legend([l2, l3], ["residuals", r"$\mathcal{K}_l$"], loc="upper right", ncol=2, prop={'size':10})
fig.text(0.03, 0.7, r"$f_\lambda \quad [\textrm{erg}\;\textrm{cm}^{-2}\;\textrm{s}^{-1}\;\textrm{\AA}^{-1}]$", rotation="vertical")
fig.subplots_adjust(hspace=0, left=0.12, right=0.95, bottom=0.17, top=0.93)
fig.savefig("../../plots/regions_23.png")
fig.savefig("../../plots/regions_23.pdf")
fig.savefig("../../plots/regions_23.svg")
plt.show()
In [ ]: